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Holographic optical traps use the forces exerted by computer-generated holograms to trap, move 
and otherwise transform mesoscopically textured materials. This article introduces methods for 
optimizing holographic optical traps' efficiency and accuracy, and an optimal statistical approach for 
characterizing their performance. This combination makes possible real-time adaptive optimization. 



A single laser beam brought to a focus with a strongly 
converging lens forms a type of optical trap widely known 
as an optical tweezer pj . Multiple beams of light passing 
simultaneously through the lens' input pupil yield mul- 
tiple optical tweezers, each at a location determined by 
its beam's angle of incidence and degree of collimation at 
the input pupil. The trap-forming laser beams form an 
interference pattern as they pass through the input pupil, 
whose amplitude and phase corrugations characterize the 
downstream trapping pattern. Imposing the same mod- 
ulations on a single incident beam at the input pupil 
would yield the same pattern of traps. Such wavefront 
modification can be performed by a computer-designed 
diffractive optical element (DOE) commonly known as a 
hologram. 

Holographic optical trapping (HOT) uses methods of 
computer-generated holography (CGH) toproject arbi- 
trary configurations of optical traps El 0, 0, El , with 
manifold applications in the physical and biological sci- 
ences as well as in industry y\. This flexible approach 
to manipulating and transforming mesoscopic matter has 
been used to assemble two- and three-dimensional struc- 
tures, to sort objects ranging in size from nanoclusters to 
living cells, and to create all-optical microfluidic pumps 
and mixers. 

This article describes refinements of the HOT tech- 
nique that help to optimize the traps' performance, as 
well as a statistically optimal analysis for rapidly char- 
acterizing them. Section [fl describes modifications to the 
basic HOT optical train that minimize defects due to 
limitations of practical implementations. Section [H] dis- 
cusses a direct search algorithm for HOT DOE computa- 
tion that is both faster and more accurate than iterative 
refinement algorithms. A method for rapidly character- 
izing each trap in a holographic array is presented in 
Section IIIII The optimal statistical methods on which 
this characterization technique is based lends itself natu- 
rally to digital video analysis of optically trapped spheres 
and can be exploited for real-time optimization. Such 
adaptive optimization is demonstrated experimentally in 
Section ITVI 
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I. IMPROVED OPTICAL TRAIN 
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FIG. 1: Simplified schematic of a holographic optical tweezer 
optical train before and after modification, (a) A collimated 
beam is split into multiple beams by the DOE, each of which 
is shown here as being collimated. The diffracted beams pass 
through the input pupil of an objective lens and are focused 
into optical traps in the objective's focal plane. The un- 
diffracted portion of the beam, shown here with the dark- 
est shading, also focuses into the focal plane, (b) The input 
beam is converging as it passes through the DOE. The DOE 
collimates the diffracted beams, so that they focus into the 
focal plane, as in (a). The undiffracted beam comes to a focus 
within the coverslip bounding the sample, (c) A beam block 
can eliminate the undiffracted beam without substantially de- 
grading the optical traps. 

Figure [Ha) shows a simplified schematic of a conven- 
tional HOT implementation. A collimated laser beam is 
imprinted with a CGH and thereafter propagates as a su- 
perposition of independent beams, each with individually 
specified wavefront characteristics 0, 0|- These beams 
are relayed to the input aperture of a high-numerical- 
aperture lens, typically a microscope objective, that fo- 
cuses them into optical traps. This figure indicates a 
transmissive DOE, although comparable results are ob- 
tained with reflective DOE's. The same objective lens 
used to form the optical traps also can be used to cre- 
ate images of trapped objects. The associated illumina- 
tion and image- forming optics are omitted from Fig. 
for clarity. Practical holograms only diffract a portion 
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of the incident light into intended modes and directions. 
Some of the incident beam is not diffracted at all, and the 
undiffracted portion typically forms an unwanted trap in 
the middle of the field of view This "central spot" 
has been removed in previous studies by spatially filtering 
the diffracted beam Practical DOE's also tend to 

project spurious "ghost" traps into symmetry-dictated 
positions within the sample. Spatially filtering a large 
number of ghost traps generally is not practical, particu- 
larly in the case of dynamic holographic optical tweezers 
whose traps move freely in three dimensions. Project- 
ing holographic traps in the off-axis Fresnel geometry, 
rather than the Fraunhofer geometry described here, au- 
tomatically eliminates the central spot |f(l |. but limits 
the number of traps that can be projected, and also does 
not address the formation of ghost traps. 

Figure H[b) shows one simple improvement to the ba- 
sic HOT design that minimizes the central spot's influ- 
ence and effectively eliminates ghost traps. Here, the 
source laser beam is converging as it passes through the 
DOE. As a result, the undiffracted central spot focuses 
upstream of the objective's normal focal plane. The de- 
gree of collimation of each diffracted beam, and thus the 
axial position of the resulting trap, can be adjusted by 
incorporating wavefront-shaping phase functions into the 
hologram's design Q , thereby returning the traps to the 
focal volume. This simple expedient allows the central 
spot to be projected into the coverslip bounding a sam- 
ple, rather than the sample itself, thereby ensuring that 
the undiffracted beam lacks both the intensity and the 
gradients needed to influence a sample's dynamics. 

An additional consequence of the traps' displacement 
relative to the converging beam's focal point is that the 
majority of ghost traps are projected to the far side of this 
point, and therefore out of the sample volume altogether. 
This is a substantial improvement for processes such as 
optical fractionation [II], U^, which rely on a precisely 
specified optical potential energy landscape. 

Even though the undiffracted beam may not create an 
actual trap in this modified optical train, it still can exert 
radiation pressure on the region of the sample near the 
center of the field of view. This is a particular problem 
for large arrays of optical traps in that the central spot 
can be brighter than the intended traps. Illuminating the 
DOE with a diverging beam 0] reduces the undiffracted 
beam's influence by projecting some of its light out of 
the optical train. In a thick sample, however, this has 
the deleterious effect of projecting both the weakened 
central spot and the undiminished ghost traps into the 
sample. 

These problems all can be mitigated by placing a beam 
block as shown in Fig. ^c) in the intermediate focal 
plane within the relay optics to spatially filter the un- 
diffracted portion of the beam. Because the trap-forming 
beams come to a focus in a different plane, they are only 
slightly occluded by the beam block, even if they pass di- 
rectly along the optical axis. The effect of this occlusion 
is minimal for conventional optical tweezers and can be 



compensated by increasing their relative brightness. 

II. ITERATIVE AND DIRECT SEARCH 
ALGORITHMS FOR HOT CALCULATION 

Holographic optical tweezers' efficacy is determined by 
the quality of the trap-forming DOE, which in turn re- 
flects the performance of the algorithms used in their 
computation. Previous studies have applied holograms 
calculated by simple linear superposition of the input 
fields 0, with best results being obtained with ran- 
dom relative phases 0, @, or with variations 0, 0, @ 
on the classic Gerchberg-Saxton and Adaptive- Additive 
algorithms [lij. Despite their general efficacy, these al- 
gorithms yield traps whose relative intensities can differ 
greatly from their design values, and typically project an 
unacceptably large fraction of the input power into ghost 
traps. These problems can become acute for complicated 
three-dimensional trapping patterns, particularly when 
the same hologram also is used as a mode converter to 
project multifunctional arrays of optical traps 4,6]. This 
section describes faster and more effective algorithms for 
HOT DOE calculation based on direct search and simu- 
lated annealing. 

The holograms used for holographic optical trapping 
typically operate only on the phase of the incident beam, 
and not its amplitude. Such phase-only holograms, 
also known as kinoforms, are far more efficient than 
amplitude- modulating holograms, which necessarily di- 
vert light away from the traps. Quite general trapping 
patterns can be achieved with kinoforms because optical 
tweezers rely for their operation on intensity gradients 
but not on phase variations. The challenge is to find a 
phase pattern in the input plane that encodes the desired 
intensity pattern in the focal volume. 

According to scalar diffraction theory, the complex 
field E(r, ip), in the focal plane of a lens of focal length / 
is related to the field, u(p) exp(i(^(p)), in its input plane 
by a Fraunhofer transform, 

E(f) = J u(p) eMMp)) cxp d 2 P , (I) 

where u(p) and <p(p) are the real-valued amplitude and 
phase, respectively, of the field at position p in the in- 
put pupil, and k = 2ir/\ is the wavenumber of light of 
wavelength A. 

If u(p) is the amplitude profile of the input laser beam, 
then (f(p) is the kinoform encoding the pattern. Most 
practical DOEs, including those projected with SLMs, 
consist of an array pj of discrete phase pixels, each of 
which can impose any of P possible discrete phase shifts 
(pj G {0, . . . , 4>p-i}. The field in the focal plane due to 
such an TV-pixel DOE is, therefore, 

N 

^(r) = X)% ex P(^i)^( r 1> ( 2 ) 
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where the transfer matrix describing the light's propaga- 
tion from input plane to output plane is 

T j (f)=e X p(-i^-y (3) 

Unlike more general holograms, the desired field in the 
output plane of a holographic optical trapping system 
consists of M discrete bright spots located at r m : 

M 

E(r) = E ^), with (4) 

m— 1 

E m (r) = a m S(r- r m ) exp(i£ m ), (5) 

where a m is the relative amplitude of the m-th trap, nor- 
malized by X)m=i \ a m\ 2 = 1, and £, m is its (arbitrary) 
phase. Here 5(f) represents the amplitude profile of the 
focused beam of light, which may be treated at least 
heuristically as a two-dimensional Dirac delta function. 
The design challenge is to solve Eqs. J3J and JIJ for 
the set of phase shifts ipj yielding the desired amplitudes 
a m at the correct locations f m , given Uj and Tj(f m ). 

The Gerchberg-Saxton algorithm and its generaliza- 
tions, such as the adaptive-additive algorithm, iteratively 
solve both the forward transform described by Eqs. © 
and 10, and also its inverse, taking care at each step to 
converge the calculated amplitudes at the output plane to 
the design amplitudes and to replace the back-projected 
amplitudes, Uj, at the input plane with the laser's actual 
amplitude profile. Appropriately updating the calculated 
input and output amplitudes at each cycle can cause the 
DOE phase tfj to converge to an approximation to the 
ideal kinoform, with monotonic convergence possible for 
some variants |14| . The forward and inverse transforms 
mapping the input and output planes to each other typ- 
ically are performed by fast Fourier transform (FFT). 
Consequently, the output positions f are discretized in 
units of the Nyquist spatial frequency. The output field, 
furthermore, is calculated not only at the intended posi- 
tions of the traps, but also at the spaces between them. 
This is useful because the process not only maximizes 
the fraction of the input light diffracted into the desired 
locations, but also minimizes the intensity of stray light 
elsewhere. The down side to this is the additional com- 
putational cost for computing fields in otherwise uninter- 
esting places. 

FFT-based iterative algorithms have additional draw- 
backs for computing three-dimensional arrays of optical 
tweezers, or mixtures of more general types of traps. To 
see this, we review how a beam-splitting DOE can be gen- 
eralized to include wavefront-shaping capabilities. A di- 
verging or converging beam at the input aperture comes 
to a focus and forms a trap downstream or upstream of 
the focal plane, respectively. Its wavefront at the input 
plane is characterized by the parabolic phase profile 



where z is the focal spot's displacement along the optical 
axis relative to the lens' focal plane. This phase profile 
can be used to move an optical trap relative to the focal 
plane even if the input beam is collimated by appropri- 
ately augmenting the transfer matrix: 

T*(f)=T 3 (f)K^ (7) 

where the displacement kernel is 

K] = exp(iv>, (/»},*)). (8) 

The result, Tj(f), replaces Tj(r) as the kernel of Eq. J3J). 

Similarly, a conventional TEMqo beam can be con- 
verted into a helical mode through 

<pW = ib, (9) 

where is the azimuthal angle around the optical axis 
and i is known as the topological charge. Such corkscrew- 
like beams focus to ring-like optical traps known as op- 
tical vortices, which can exert torques as well as forces. 
The topology-transforming kernel K | = exp(iipi(pj)) can 
be composed with the transfer matrix in the same man- 
ner as the displacement-inducing K*. A variety of anal- 
ogous phase-based mode transformations have been de- 
scribed, each with applications to single-beam optical 
trapping 0. All can be implemented with appropriate 
transformation kernels. Moreover, different transforma- 
tion operations can be applied to each beam in a holo- 
graphic trapping pattern independently, resulting in gen- 
eral three-dimensional configurations of diverse types of 
optical traps |6(- 

Calculating the phase pattern ipj encoding multifunc- 
tional three-dimensional optical trapping patterns re- 
quires only a slight elaboration. The primary require- 
ment is to measure the actual intensity projected by 
ipj into the m-th trap at its focus. If the associated 
diffraction- generated beam has a non-trivial wavefront, 
then it need not create a bright spot at its focal point. On 
the other hand, if we assume that tpj creates the required 
type of beam for the m-th trap through a phase modu- 
lation described by the transformation kernel Kj^ m , then 
applying the inverse operator, Kj m , in Eq. J2J) would re- 
store the focal spot. This principle was first applied to 
creating three dimensional trap arrays Q in which sepa- 
rate translation kernels were used to project each desired 
optical tweezer back to the focal plane as an intermediate 
step in each iterative refinement cycle. Computing the 
light projected into each plane of traps in this manner 
involves a separate Fourier transform for the entire plane 

3. 

Much of this effort is eliminated by computing the field 
only at the traps' positions, as 

N 

E m (f m ) = K iX T ^m) e Mi<Pj), (10) 



(6) 



and comparing the resulting amplitudes a m = \E m \ with 
the design values @. Unlike the FFT-based approach, 
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this per-trap algorithm does not directly optimize the 
field between the traps. However, it also eliminates the 
need to account for interplane propagation. Moreover, if 
the values of a m match the design values, no light is left 
over to create ghost traps. 

The per-trap calculation suffers from its own shortcom- 
ings. The only adjustable parameters in Eqs. (JSJ) and 
(|TU)l are the relative phases £ m of the projected traps. 
These M — 1 real-valued parameters must be adjusted 
to optimize the choice of discrete- valued phase shifts, ipj , 
subject to the constraint that the amplitude profile Uj 
matches the input laser's. With so few free parameters, 
however, finding an optimal hologram is not likely. 

Equation QlOfl suggests an alternative approach for 
computing DOE functions for discrete HOT patterns. 
The operator K~^ ll Tj(r m ) describes how light in the 
mode of the m-th trap propagates from position pj on 
the DOE to the trap's projected position r m in the lens' 
focal plane. If we were to change the DOE's phase (fij 
at that point, then the superposition of rays composing 
the field at r m would be affected. If the change improves 
the overall trapping pattern, then we would be inclined 
to retain it, and seek other such improvements. If, in- 
stead, the result were less good, we would restore (pj to 
its former value and look elsewhere. This is the basis for 
direct search algorithms. 

The simplest direct search involves selecting a pixel at 
random from a trial phase pattern, changing its value to 
any of the P — 1 alternatives, and computing the effect on 
the projected field. This operation can be performed ef- 
ficiently by calculating only the changes at the M traps' 
positions. The updated trial amplitudes then are com- 
pared with their design values and the proposed change 
is accepted if the overall error is reduced. The process is 
repeated until the result converges to the design or the 
acceptance rate for proposed changes dwindles. 

A successful and efficient direct search requires an ef- 
fective assessment of errors. The standard cost func- 
tion, x 2 — Em=i(4 _ e^Tjf^) 2 , assesses the mean- 
squared deviations of the m-th trap's projected inten- 
sity I rn — |a m | 2 from its design value Im \ assuming an 
overall diffraction efficiency of e. It requires an accurate 
estimate for e and places no emphasis on uniformity in 
the projected traps' intensities. An alternative proposed 
by Meister and Winfield 

C = -(I) + fa, (11) 

avoids both shortcomings. Here, (I) is the mean intensity 
at the traps and 
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(12) 



measures the deviation from uniform convergence to the 
design intensities. Selecting 
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minimizes the total error and accounts for non-ideal 
diffraction efficiency. The weighting fraction / sets the 
relative importance attached to diffraction efficiency ver- 
sus uniformity. 

In the simplest direct search for an optimal phase dis- 
tribution, any candidate change that reduces C is ac- 
cepted, and all others are rejected. In a worst-case im- 
plementation, the number of trials required for practi- 
cal convergence should scale as NP, the product of the 
number of phase pixels and the number of possible phase 
values. In practice, this rough estimate is accurate if P 
and N are comparatively small and if the starting phase 
function is either uniform or purely random. Much faster 
convergence can be obtained by starting from the phase 
obtained by superposing the desired beams with random 
relative phases and ignoring amplitude modulations. In 
this case, convergence typically is obtained within N tri- 
als, even for fairly complex trapping patterns, and thus 
requires a computational effort comparable to the initial 
superposition. 

As a practical demonstration, we implement a chal- 
lenging quasiperiodic array of optical traps. The traps 
are focused with a 100 x NA 1.4 S-Plan Apo oil immersion 
objective lens mounted in a Nikon TE-2000U inverted 
optical microscope. The traps are powered by a Co- 
herent Verdi frequency-doubled diode-pumped solid state 
laser operating at a wavelength of 532 ran. Computer- 
generated phase holograms are imprinted on the beam 
with a Hamamatsu X8267-16 parallel-aligned nematic 
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FIG. 2: (a) Design for 119 identical optical traps in a two- 
dimensional quasiperiodic array, (b) Trapping pattern pro- 
jected without optimizations using the adaptive-additive al- 
gorithm, (c) Trapping pattern projected with optimized op- 
tics and adaptively corrected direct search algorithm, (d) 
Bright-field image of colloidal silica spheres 1.53 fim in di- 
ameter dispersed in water and organized in the optical trap 
array. The scale bar indicates 10 /an 
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liquid crystal spatial light modulator (SLM). This SLM 
can impose phase shifts up to 2tt radians at each pixel in 
a 768 x 768 array. The face of the SLM is imaged onto the 
objective's 5 mm diameter input pupil using relay optics 
designed to minimize aberrations. The beam is directed 
into the objective with a dichroic beamsplitter (Chroma 
Technologies), which allows images to pass through to a 
low-noise charge-coupled device (CCD) camera (NEC TI- 
324AII). The video stream is recorded as uncompressed 
digital video with a Pioneer 520H digital video recorder 
(DVR) for processing. 

Figure Ufa) shows the intended planar arrangement 
of 119 holographic optical traps. Even after adaptive- 
additive refinement, the hologram resulting from simple 
superposition with random phase fares poorly for this 
aperiodic pattern. Figure[2Ib) shows the intensity of light 
reflected by a front-surface mirror placed in the sample 
plane. This image reveals extraneous ghost traps, an ex- 
ceptionally bright central spot, and large variability in 
the intended traps' intensities. Imaging photometry on 
this and equivalent images produced with different ran- 
dom relative phases for the beams yields a typical root- 
mean-square (RMS) variation of more than 50 percent in 
the projected traps' brightness. The image in Fig. |2fc) 
was produced using the modified optical train described 
in Sec. U and the direct search algorithm described in 
Sec. [n] and suffers from none of these defects. Both the 
ghost traps and the central spot are suppressed, and the 
apparent relative brightness variations are smaller than 5 
percent, a factor often improvement. Figure |2Id) shows 
119 colloidal silica spheres, 2a — 1.5±0.3 /zm in diameter 
(Bangs Labs, lot 5238), dispersed in water at T = 27°C 
and trapped in the quasiperiodic array. 

To place the benefits of the direct search algorithm on 
a more quantitative basis, we augment standard figures 
of merit with those introduced in Ref. [l^. In particular, 
the DOE's theoretical diffraction efficiency is commonly 
defined as 
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and its root-mean-square (RMS) error as 

a 



max(7 m ) 



(14) 



(15) 



The resulting pattern's departure from uniformity is use- 
fully gauged as [l5| 



max(/ m /7, 



max(J m + min(7 m /I, 



(16) 



Figure |3 shows results for a HOT DOE encoding 51 
traps, including 12 optical vortices of topological charge 
I = 8, arrayed in three planes relative to the focal plane. 
The excellent results in Fig. were obtained with a sin- 
gle pass of direct-search refinement. The resulting traps, 
shown in the bottom three images, again vary from their 



planned relative intensities by less than 5 percent. In this 
case, the spatially extended vortices were made as bright 
as the point-like optical tweezers by increasing their re- 
quested relative brightness by a factor of 15. This single 
hologram, therefore, demonstrates independent control 
over three-dimensional position, wavefront topology, and 
brightness of all the traps. Performance metrics for the 
calculation are plotted in Fig. |3Ib) as a function of the 
number of accepted single-pixel changes, with an overall 
acceptance rate of 16 percent. Direct search refinement 
achieves greatly improved fidelity to design over ran- 
domly phase superposition at the cost of a small fraction 
of the diffraction efficiency and roughly doubled compu- 
tation time. The entire calculation can be completed in 
the refresh interval of a typical liquid crystal spatial light 
modulator. 



III. OPTIMAL CHARACTERIZATION 

At least numerically, direct search algorithms are both 
faster and better at calculating trap-forming DOEs than 
iterative refinement algorithms. The real test, however, 
is in the projected traps' ability to trap particles. A 
variety of approaches have been developed for gauging 
the forces exerted by optical traps. The earliest involved 
measuring the hydrodynamic drag required to dislodge 
a trapped particle 16]. This has several disadvantages, 
most notably that it identifies only the marginal escape 
force in a given direction and not the trap's actual three- 
dimensional potential. Most implementations, further- 
more, failed to collect sufficient statistics to account for 
thermal fluctuations' role in the escape process, and did 
not account adequately for hydrodynamic coupling to 
bounding surfaces. 

Complementary information can be obtained by mea- 
suring a particle's thermally driven motions in the trap's 
potential well 0,0, El- For instance, the measured 
probability density P(r) for displacements r is related to 
the trap's potential V(r) through the Boltzmann distri- 
bution 



P(r) cx exp(-/?V(r)), 



(17) 



where /3 _1 = ksT is the thermal energy scale at tem- 
perature T. Similarly, the power spectrum of f(t) for a 
harmonically bound particle is a Lorentzian whose width 
is the viscous relaxation time of the particle in the well 

dm. 

Both of these approaches require amassing enough 
data to characterize the trapped particle's least proba- 
ble displacements, and therefore oversample the trajec- 
tories. Oversampling is acceptable when data from a 
single optical trap can be collected rapidly, for example 
with a quadrant photodiode [13, EE 111. E3 • Tracking 
multiple particles in holographic optical traps, however, 
is most readily accomplished through digital video mi- 
croscopy [22J, which yields data much more slowly. An- 
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FIG. 3: A three-dimensional multifunctional holographic optical trap array created with the direct search algorithm, (a) 
Refined DOE phase pattern, (b), (c) and (d) The projected optical trap array at z = —10 /im, /im and +10 /jm. Traps are 
spaced by 1.2 fim in the plane, and the 12 traps in the middle plane consist of I = 8 optical vortices, (e) Performance metrics 
for the hologram in (a) as a function of the number of accepted single-pixel changes. Data include the DOE's overall diffraction 
efficiency as defined by Eq. 1141 . the projected pattern's RMS error from Eq. (1151 . and its uniformity, 1 — u, where u is defined 
in Eq. (17)1 . 



alyzing video data with optimal statistics [23| offers the 
benefits of thermal calibration by avoiding oversampling. 

An optical trap is accurately modeled as a harmonic 
potential energy well [3 El Efl E3 , 



V(f) 



1 V 



(18) 



i=l 



with a different characteristic curvature m along each 
axis. This separable form admits a one-dimensional anal- 
ysis. The trajectory of a colloidal particle localized in a 
viscous fluid by a harmonic well is described by the one- 
dimensional Langevin equation [24[ 



T 



(19) 



where the autocorrelation time r = j/k, is set by the 
viscous drag coefficient 7 and the curvature of the well, 
k, and where £(t) describes Gaussian random thermal 
forcing with zero mean, (£(£)) = 0, and variance 



(£(*)£(*)> = ^ *(*-«)• 

7 



(20) 



If the particle is at position xq at time t = 0, its trajectory 
at later times is given by 



x(t) — xq exp 



t 



£(s) exp 



t - s 



ds. (21) 



Sampling such a trajectory at discrete times tj = j At, 
yields 



Xj+i = xq exp 



t 



£(s) exp - 



tj+l - s 



ds 



+ / £(a) exp - 



tj+i - s 



ds 



(22) 



At 



= 4>Xj + aj+i, where <fi = exp I J , (23) 

and where Oj+i is a Gaussian random variable with zero 
mean and variance 



1 — exp 



(24) 



Because <f> < 1, Eq. (|23|l is an example of an autoregres- 
sive process |23|, which is readily invertible. In principle, 
the particle's trajectory {xj} can be analyzed to extract 
4> and cr^, and thus the trap's stiffness, K, and the parti- 
cle's viscous drag coefficient 7. 

In practice, however, the experimentally measured par- 
ticle positions yj differ from the actual positions Xj by 
random errors bj, which we assume to be taken from a 
Gaussian distribution with zero mean and variance a\. 
The measurement then is described by the coupled equa- 
tions 



<j)Xj 



and 



y.i 



(25) 
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where bj is independent of a 3 . We still can estimate <f> and 
a\ from a set of measurements {yj} by first constructing 
the joint probability 



(N + 1) matrix 
/1 + 2 



JY 



p{{xi},{ yi }\(j),al,al) = Y[ 
3=2 
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N 
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(26) 
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2^ 
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Physically, this involves imparting an impulse, (Zjv+ij 
that translates the particle from its last position, xn, 
to its first, x\. Because diffusion in a potential well is a 
stationary process, the effect of this change is inversely 
proportional to the number of measurements, N. 

With this approximation, the determinant and inverse 
of Ad, are given by 



(27) 



N r i 

det^ = n - 



The probability density for measuring the trajectory 
{yj}, is then the marginal 



1 



— 2(f> cos 



(2im\ 



V N J 



(33) 



and 



iV 



P{{yj}\<t>,°-l,°b) = \ P{{xj},{y 3 }\&°l,Vb)dx 1 ---dx N {A ctj 1 ) a (3 = Jr 
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so that the conditional probability for the measured tra- 
jectory, {yj}, is 



4^ 



p{{y 3 }\<P,alal) = {2^y 



(29) 
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where y = (j/i, . . . , yjv), (y) T is its transpose, I is the 
N x N identity matrix, and 
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with the tridiagonal memory tensor 
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where y m is the m-th component of the discrete Fourier 
transform of {y n }- This can be inverted to obtain the 
likelihood function for cj>, of, and of: 



(31) 



,L^, ( r 2 ,a 2 |{y s ;}) = -y ln27r 
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Calculating the determinant, det(A^), and inverse, A^ 1 , 
of A<p is greatly facilitated if we artificially impose time 
translation invariance by replacing with the (N+l) x 



71=1 



— 20 cos 



27T7J 



(36) 



Best estimates ((f), cr%,cr%) for the parameters (</>, of, of) 
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are solutions of the coupled equations 
dL dL dL 



dal dcr? 



= 0. 



(37) 



A. Case 1: No measurement errors (of = 0) 

Equations (|37|) can be solved in closed form if cr^ = 0. 
In this case, 



Centroid resolutions of 10 nm or better can be attained 
routinely for micrometer-scale colloidal particles in wa- 
ter using conventional bright-field imaging. In practice, 
however, mechanical vibrations, video jitter and other 
processes may increase the measurement error. Quite 
often, the overall measurement error is most easily as- 
sessed by increasing the laser power to the optical traps 
to minimize the particles' thermally driven motions. In 
this case, yj « bj 



and of can be estimated directly. 



lo= — and al = c 
Co 



where 



1-1^ 

Co 



(38) 



C. Trap characterization 

The stiffness and viscous drag coefficient can be esti- 
mated simultaneously as 



1 N 



m) mod N 



(39) 



is the barrel autocorrelation of {yj} at lag m. The asso- 
ciated statistical uncertainties are 



A0 O 



Nc > 



and Acr^n = a% n ■ 



(40) 



In the absence of measurement errors, Co and Ci consti- 
tute sufficient statistics for the time series |23| and thus 
embody all of the information that can be extracted. 



B. Case 2: Small measurement errors (of -C of) 

The analysis is less straightforward when of ^ be- 
cause Eqs. (|37f) no longer are simply separable. The sys- 
tem of equations can be solved approximately if of <C of . 
In this case, the best estimates for the parameters can be 
expressed in terms of the error-free estimates as 
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'o0 



co 



00 — 

CO 



(41) 



to first order in <j\I<j\, with statistical uncertainties prop- 
agated in the conventional manner. Expansions to higher 
order in a\ /erf involve additional correlations, and the 
exact solution involves correlations at all lags m. If mea- 
surement errors are small enough for Eq. I|41|) to apply, 
the computational savings relative to other approaches 
can be substantial, and the amount of data required to 
achieve a desired level of accuracy in the physically rele- 
vant quantities, k and 7, can be reduced dramatically. 

The errors in locating colloidal particles' centroids can 
be calculated from knowledge of the images' signal to 
noise ratio and the optical train's magnification |22j . 



K 1 — 



k R T 



and 



7 



1 - 



b± at 



k B TAt a 2 in 



(42) 



with error estimates, An and A7, given by 
2 




If the measurement interval, At, is much longer than the 
viscous relaxation time r = 7/K, then e/> vanishes and 
the error in the drag coefficient diverges. Conversely, 
if At is much smaller than r, then <j) approaches unity 
and both errors diverge. Consequently, this approach 
does not benefit from excessively fast sampling. Rather, 
it relies on accurate particle tracking to minimize A(j) 
and Acr^. For trap-particle combinations with viscous 
relaxation times exceeding a few milliseconds and typ- 
ical particle excursions of at least 10 nm, digital video 
microscopy provides the resolution needed to simultane- 
ously characterize multiple optical traps 22]. 

In the event that measurement errors can be ignored 
(of < 



up 

and 

7o 
k B TAt 



1 

co 




co In fa 



1 ± 



A70 
7o 



(45) 



where 



N 



A70 

7o 



2cf ln(a)-c 



(46) 
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FIG. 4: Power dependence of (a) the trap stiffness, (b) the viscous drag coefficient and (c) the viscous relaxation time for a 
1.53 /xm diameter silica sphere trapped by an optical tweezer in water. 



These results are not reliable if c\ < of, which arises 
when the sampling interval At is much longer or much 
shorter than the viscous relaxation time, r. Accurate 
estimates for n and 7 still may be obtained in this case 
by applying Eq. 

As a practical demonstration of this formalism, we 
analyzed the thermally driven motions of a single sil- 
ica sphere of diameter 1.53 /jm (Bangs Labs lot number 
5328) dispersed in water and trapped in a conventional 
optical tweezer. With the trajectory resolved to within 
6 nm at 1/30 sec intervals, 1 minute of data suffices to 
measure both k and 7 to within 1 percent error using 
Eqs. 142fl . The results plotted in Fig. Ufa) indicate trap- 
ping efficiencies of k x /P — K y /P — 142 ± 3 pN//imW. 
Unlike k, which depends principally on c , 7 also depends 
on ci, which is accurately measured only for r > 1. Over 
the range of laser powers for which this condition holds, 
we obtain the expected 7^/70 = ly/la = 1-0 ± 0.1, as 
shown in Fig. EJb) ■ The viscous relaxation time becomes 
substantially shorter than our sampling time for higher 
powers, so that estimates for 7 and for the error in 7 
both become unreliable, as expected. 



IV. ADAPTIVE OPTIMIZATION 

Optimal statistical analysis offers insights not only into 
the traps' properties, but also into the properties of the 
trapped particles and the surrounding medium. For ex- 
ample, if a spherical probe particle is immersed in a 
medium of viscosity 77 far from any bounding surfaces, 
its hydrodynamic radius a can be assessed from the mea- 
sured drag coefficient using the Stokes result 7 = 6^770. 
The viscous drag coefficients, moreover, provide insights 
into the particles' coupling to each other and to their 
environment. The independently assessed values of the 
traps' stiffnesses then can serve as a self-calibration in 



microrheological measurements and in measurements of 
colloidal many-body hydrodynamic coupling [25j. In 
cases where the traps themselves must be calibrated ac- 
curately, knowledge of the probe particles' differing prop- 
erties gauged from measurements of 7 can be used to dis- 
tinguish variations in the traps' intrinsic properties from 
variations due to differences among the probe particles. 

These measurements, moreover, can be performed 
rapidly enough, even at conventional video sampling 
rates, to permit real-time adaptive optimization of the 
traps' properties. Each trap's stiffness is roughly propor- 
tional to its brightness. So, if the m-th trap in an array 
is intended to receive a fraction |a m | 2 of the projected 
light, then instrumental deviations can be corrected by 
recalculating the CGH with modified amplitudes: 



a m ^a m \ ^ l=1 - . (47) 

Analogous results can be derived for optimization on the 
basis of other performance metrics. A quasiperiodic pat- 
tern similar to that in Fig.^c) was adaptively optimized 
for uniform brightness, with a single optimization cycle 
yielding better than 12 percent variance from the mean. 
Applying Eqs. to data from images such as Fig.^d) 
allows us to correlate the traps' appearance to their ac- 
tual performance. 

With each trap powered by 3.4 mW, the mean viscous 
relaxation time is found to be r/ At = 1.14 ± 0.11. We 
expect reliable estimates for the viscous drag coefficient 
under these conditions, and the result 7/70 = 0.95 ±0.10 
with an overall measurement error of 0.01, is consistent 
with the manufacturer's rated 10 percent polydispersity 
in particle radius. Variations in the measured stiffnesses, 
(k x ) = 0.38±0.06 pN//xm and («„) = 0.35±0.10 pN//xm, 
can be ascribed to a combination of the particles' poly- 
dispersity and the traps' inherent brightness variations. 
This demonstrates that adaptive optimization based on 
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the traps' measured intensities also optimizes their per- 
formance in trapping particles. 

V. SUMMARY 

The quality and uniformity of the holographic opti- 
cal traps projected with the methods described in the 
previous sections represent a substantial advance over 
previously reported results. We have demonstrated that 
optimized and adaptively optimized HOT arrays can be 
used to craft highly structured potential energy land- 
scapes with excellent fidelity to design. These optimized 



landscapes have potentially wide-ranging applications in 
sorting mesoscopic fluid-borne objects through optical 
fractionation 0, in fundamental studies of trans- 
port [U |2(|, dynamics [53, HU and phase transitions 
in macromolecular systems, and also in precision holo- 
graphic manufacturing. 
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